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Abstract 

ON 

Solving large polynomial systems with coefficient parameters are ubiquitous and 
constitute an important class of problems. We demonstrate the computational power 
of two methods-a symbolic one called the Comprehensive Grobner basis and a nu- 
merical one called the cheater's homotopy-applied to studying both potential energy 
landscapes and a variety of questions arising from geometry and phenomenology. 
Particular attention is paid to an example in flux compactification where important 
physical quantities such as the gravitino and moduli masses and the string coupling 
can be efficiently extracted. 
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1 Introduction 

The hypersurface defined by a potential energy function V(x; a), with X — • • • 5 % iv) 
being the variables and a — (ax, ... , at) being a set of parameters, is called the poten- 
tial energy landscape (PEL) of a given physical model. The landscape, in particular, 
refers to the potentially large number of parameters which constitutes a moduli space. 
The special points of a PEL, defined as the critical (or stationary) points of = 0, 
with 1 < i < N, give crucial information about the physical system depending on 
the problem at hand. 

In recent years, a huge influx of advancements have come through from many 
areas in theoretical physics and chemistry in understanding the PEL and its rela- 
tion to various physical and chemical properties. The research areas where the PEL 
methods have have been immensely successful in explaining the underlying physics or 
chemistry include clusters [1-6], disordered systems and glasses [7, 8], biomolecules, 
protein folding, string phenomenology [9], and within this area flux compactifica- 
tions [10-16]. 

Due to the importance of these critical points of the PELs, finding them in real- 
istic models has been an active research area for quite some time. The stationary 
equations for any realistic model are invariably non-linear which are known to be 
extremely difficult to solve in general. Various numerical techniques exist based on 
the Newton-Raphson method and its sophisticated variants [17-21] where a random 
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initial guess is refined to attain a single solution of the system. However, in all these 
methods, even after feeding a large number of random initial guesses, one can never 
be sure of obtaining all the solutions. 

While many solutions as opposed to all solutions may be fine in many applications, 
often one has to find all the solutions. The problem of not obtaining all the solutions 
becomes crucial when one needs the information about all the local minima or the 
global minimum, etc. For example, in string theory models, where in many cases the 
potential energy landscape is defined by an effective four-dimensional, N = 1 super- 
gravity scalar potential V(K, W) given a Kahler potential K, and a superpotential 
W, one has to find all the minima. These minima are called string vacua and ad- 
dressing the plethora of solutions is one of the most important of current theoretical 
challenges. 

Fortunately, in a large class of models, the equations are polynomials, or at least 
they have polynomial-like non-linearity. 1 Once they are identified as a system of 
polynomial equations, we can use techniques from algebraic geometry, or more specif- 
ically, computational algebraic geometry to solve them by systematically transform- 
ing a given system of multivariate polynomial equations to another system which is 
easier to solve and has the same solutions as the original one. The new system of 
equations is called a Grobner basis (GB), and the algorithm to compute one is called 
the Buchberger algorithm (BA). 

The biggest advantage here is that one can find all the solutions of the given 
system once the computation is finished. Computational algebraic geometry, which 
is essentially a set of techniques based on the Grobner basis method, has become 
one of the most useful tools to study a number of phenomena in theoretical physics. 
Recently, the rich interplay between algebraic geometry and theoretical physics, espe- 
cially in gauge and string theory, has been an active area of research [23]. Activities 
in these areas have been enhanced with the increased power of computers and the 
development of algorithms in computational algebraic geometry. 

More specifically, a variety of methods have been used to study the moduli space 
of vacua over the past few years [24, 25] based on symbolic computational algebraic 
geometry, most of whose sub-methods and sub-algorithms rely on Grobner basis 
techniques (cf. [26] for an overview on the method). For convenience, a freely avail- 
able computational package, StringVacua, which is a Mathematica package specifically 
designed for phenomenologists [22, 27, 28] exists. StringVacua interfaces with the ad- 
vanced computational algebraic geometry package Singular [29]. Using StringVacua, 
one can extract important information such as the dimension of the vacuum, the 
number of real roots in the system, stability and supersymmetry of the potential, 

1 In the presence of non-perturbativc effects where exponential terms contribute, one could still 
reduce the system into polynomials of Lambert functions, perform standard polynomials manipu- 
lations and treat the Lambert functions numerically [22]. 
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or the branches of moduli space of vacua, etc., using a regular desktop machine in 
many circumstances. 

However, the GB method is known to suffer from exponential complexity, i.e., the 
computation time and the RAM required by the BA algorithm increases exponen- 
tially with the number of variables, equations, degree, and terms in each polynomial; 
it is usually less efficient for systems with irrational coefficient parameters a; and it 
is also highly sequential, i.e., very difficult to parallelize the algorithm and put on a 
big cluster. 

To overcome these shortcomings of the GB methods, a different approach called 
numerical algebraic geometry (NAG) was recently introduced. Its core algorithm, 
called numerical polynomial homotopy continuation (NPHC), guarantees to find all 
of the solutions unlike other numerical methods such as the Newton-Raphson and 
its variants. Moreover, unlike the Grobner basis techniques, the NPHC method is 
"embarrassingly parallelizable" , and hence one can solve more complicated systems 
efficiently using computer clusters. 

NAG was introduced in particle theory and statistical mechanics areas in Ref . [30] . 
Subsequently, the NPHC method was used to solve systems arising in numerous phys- 
ical phenomena in lattice field theories [31, 32], statistical physics [33-36], particle 
phenomenology [37, 38], and string phenomenology [39-41]. 

1.1 Parametric Potentials 

As mentioned above, in generic physical applications, the potential energy function is 
defined over a possibly vast parameter-space, where each point a represents a differ- 
ent physical situation. For example, in a statistical mechanics model, the parameters 
are the disorders [7]. Another example in theoretical chemistry where the models are 
Morse clusters, the parameters represent the strength of the inter-particle potential 
[42-45]. As a third example, in lattice field theory the parameters represent different 
background fields in the models [31, 32, 46, 47]. Another situation is in flux com- 
pactifications within string phenomenology where the parameters represent the flux 
quanta. 

The flux quanta are discrete parameters that are given as integrals of n-forms along 
n-cycles in a compact space, see §4. Their discreteness arises as a Dirac quantization 
condition. For many compactification manifolds, for instance Calabi-Yau spaces, 
these cycles exist on the order of 10 to 100 where the flux quanta can be chosen 
independently. This yields an exponentially large parameter space that is only limited 
by conditions of conservation of certain charges in the compact space. 

The parametric potentials adds one more, rather severe, problem where as the 
parameters enter the system of equations, one has to solve the system of parametric 
equations. This is in fact a quite difficult problem. If we use the above mentioned 
methods directly, we have to specify numerical values of the parameters before using 
the methods. In other words, we have to solve the system separately from scratch 
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for each point in the parameter space. In practice, due to the large number of phys- 
ically interesting parameter-points, this crude way of solving the problem becomes 
prohibitively time consuming as well as computationally expensive. Traditionally, 
one resorts to assigning generic, random values to these parameters. However, this 
misses the intricacies of special parameters. 

In this paper, we introduce two methods which deal with parametric systems ex- 
tremely efficiently The first method is called the Comprehensive Grobner basis 
(CGB) method. This is a completely symbolic method where given a parametric 
system of equations, CGB yields a new system, leaving all the parameters in the 
symbolic form, which is a GB for all values of the parameters and also for all the 
specializations, i.e., special values of the parameters. Once we obtain a CGB for a 
given system, it then only amounts to inputting the specific values of parameters and 
extracting the corresponding solutions out. Thus, we can efficiently solve the system 
for as many parameter-points as we desire since the system is solved for the whole 
parameter space. However, as one may surmise, this method has the same afore- 
mentioned short-comings as the usual GB method. Nevertheless, when successful in 
finishing the computation, the CGB method reduces the amount of work in finding 
string vacua over a large space of flux parameters drastically. 

The second method is called the cheater's homotopy (CH) method, which is 
based on NAG. In NPHC, one first has to estimate an upper bound on the number 
of solutions of the given system so that one can then construct another system which 
has the same number of solutions as the estimate. To reduce the computation in this 
method, we must come up with a tighter upper bound on the number of solutions. 
Three of the most important solution bounds were demonstrated in our previous 
works [30, 31, 33, 40]. However, in most of the realistic systems, the sparsity of the 
systems (i.e., the number of monomials in each equation may be only a few) is not 
fully taken into account. Thus, even though the original system may have only a few 
solutions, we may have to track many paths making it computationally expensive. 
Cheater's homotopy relies on the fact that the maximum number of solutions of a 
parametric system of polynomial equations over all the parameter-points is that at 
a generic parameter-point. The number of paths to be tracked is usually drastically 
smaller than any of the other upper bounds for the cheater's homotopy [48]. 

The second, and perhaps the most crucial, advantage is that in the cheater's 
homotopy we obtain the start system once for for a generic parametric system, i.e., 
we then can use the same start system for arbitrary parameter-points in the entire 
parameter-space. Several more advantages of this method exist when solving these 
systems on a computer cluster, but most importantly the method is "embarrassingly 
parallelizable" , a subject matter to which we shall later return. A very efficient 
package, which is yet to be publicly available, called Paramotopy [49], uses all the 
computational advantages of the method in our favor and can deal with hundreds of 
thousands of parameter-points. In this paper, we heavily rely on Paramotopy for our 



5 



computations. 

We would like to make the clarification that the above two methods are based 
on complex algebraic geometry which means that the variables are first considered 
to be complex, and then once the solutions are obtained, only purely real solutions 
are retained for analysis if the system was supposed to be made of real variables. 
However, there is a more recent method from real algebraic geometry, based on 
the discriminant varieties [50, 51], which treats both variables and parameters as 
reals from the beginning and gives a different set of information than the methods 
described in this paper. The details on this method with applications arising from 
physics will be addressed in forthcoming works. 

We first introduce the Comprehensive Grobner basis method and solve a few toy 
models in §2. Then, to make the paper self-contained, we briefly explain the NPHC 
method before explaining the cheater's homotopy method in §3. Therein, we also 
introduce the salient features of the Paramotopy package. In §4, we solve several toy 
models as well as more realistic models such as flux compactifications on the quintic 
manifold. We will extract some interesting physics using the string vacua of these 
models over a large space of fluxes. Finally, we conclude with remarks and prospects 
in §5. 

2 Algebraic Geometry &; the Comprehensive Grobner Basis 
Method 

In this section, we explain the concept of a Comprehensive Grobner basis in more 
detail and some basics of algebraic geometry to set the nomenclature and also to 
facilitate the reader. We first introduce few technical terms leading to the definition 
of a Grobner basis. Then, we will explain what a Comprehensive Grobner basis is. 
Readers uninterested in the technicalities involved in the definition may freely skip 
subsection §2.1 after reading the next two paragraphs. 

Roughly speaking, given a system of polynomial equations, the Buchberger Algo- 
rithm (BA), or its refined variants, compute a new equivalent system of polynomials, 
called a Grobner basis [52] which has nicer properties; using the BA on multivariate 
polynomial systems is analogous to Gaussian elimination for linear systems. Nowa- 
days, efficient variants of the BA are available, e.g., F4 [53], F5 [54], and Involution 
Algorithms [55]. Symbolic computation packages such as Mathematica, Maple, Re- 
duce, etc., have built-in commands to calculate a Grobner basis. Moreover, Singular 
[29], COCOA [56], and Macaulay2 [57] are specialized packages for computational al- 
gebraic geometry available freely. MAGMA [58] is also such a specialized package 
available commercially. 

A system of polynomial equations with parameters is called a parametric system; 
finding the critical points of a polynomial PEL is precisely such a system. If one 
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is interested in solving the system at finitely many points on the parameter-space, 
then inserting the numerical values of the parameters in the system and obtaining a 
Grobner basis is a quick escape, especially for a small number of parameters. A nat- 
ural question to ask is if it is possible to obtain a Grobner basis for a given monomial 
ordering in terms of the symbolic form of the parametric coefficients, valid for all 
its special cases, called specializations. One can indeed compute such a 'parametric 
Grobner basis' called Comprehensive Grobner basis (CGB) [59]. Algorithmically, we 
use the internal libraries of Singular to compute the CGB in this paper. 

2.1 The Comprehensive Grobner Basis 

The technicalities in this subsection will lead to a very useful result, namely that 
we can transform a given system of multivariate polynomial equations to another 
one which has the same solutions but is easier to solve. Here, the original system 
is considered as a basis of an algebraic object, called an ideal. Then an important 
result, that an appropriate change of this basis leaves the solution space unchanged, 
is used. 

Polynomial Rings 

We define a polynomial / as 



Here, the sum is over a finite number of m-tuples a = («i,...,a m ) and x a = 
monomial with all a, being non-negative integers. The coefficients 
a a and the variables Xi take values from the field K. However, to use the results of 
Algebraic Geometry to its full extent, unless otherwise specified, we will take K — C. 

Now, if K[ X\, . . . ,x m ] is the set of all polynomials in variables x±,..., x m with 
coefficients in K, then / e K [xi, . . . , x m ] can be viewed as a function / : K m — > K 
where K m is the affine space of all coefficients. Thus, the sum and product of 
two polynomials is a polynomial, and a polynomial / divides a polynomial g if 
and only if g = fh, for some h G K[xi, . . . ,x m }. Using this, it can be shown 
that under addition and multiplication, K[x±, . . . ,x m ] satisfies all of the field axioms 
except for the existence of a multiplicative inverse because - is not a polynomial. 
Indeed, K[xi, . . . , x m ] satisfies the axioms for a commutative ring, or more precisely 
a polynomial ring 2 . 



One can now view all the polynomials of a system of polynomial equations as elements 
of a polynomial ring. Hence, one can also define a corresponding vector space, called 
an ideal. More specifically, an ideal J is a subset of K[ 1 with the following; 

properties: 

2 For a nice discussion on related topics, the reader is referred to [52, 60]. 





Ideal 
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1. G /, 



2. / + g E I for all f,g e I, and 

3. hf E I for f E I and ft. 6 Afxi . . . , x m ]. 

Consider any f E I G K[ ]. If / can be written as / = J2 a a ah a with 

a a E K and /i Q G if [#i, • • • , x m ], then we write 7 = C if[:£i, . . . , x m ]. If the 

indexing set a is finite, say with cardinality t, then I is called a finitely generated 
ideal. The polynomials hi, ... , h t are then said to be a finite basis of i, and we write 
I = (hi, h t ). 

Affine Variety 

So far, we have introduced the algebraic counterpart of Algebraic Geometry. The 
solution space of a given ideal is called a variety. Specifically, an affine variety of an 
ideal I = (hi, . . . , h t ) is the set of common zeros of polynomials hi . . . , h t in affine 
space, denoted as V(hi, . . . , h t ) or V(I). 

Grobner Basis 

The formalism of Algebraic Geometry turns out to be very helpful. Interpreting the 
polynomials hi as a basis of /, we can change the basis to, say, (Hi, . . . , H s ). Then, 
it can be shown that the solution space remains unchanged in an appropriate change 
of basis, that is, V(hi, . . . , h t ) = V(Hi, . . . , H s ). In essence, we use computational 
techniques to find a basis that is easier to deal with than the original one, in a certain 
sense. Such a basis is called a Grobner basis. 

In linear algebra, such a change of basis can be done via Gaussian Elimination, and 
the new basis is the familiar Row-Echelon form. In general, an algorithm to obtain 
a Grobner basis performs a specific set of algebraic operations including factorizing 
and dividing the polynomials. In any algorithm that computes a Grobner basis, the 
division requires one to impose a total order among the monomials. This is called a 
monomial ordering. 

A monomial ordering is a relation V on the set of monomials x a , a E Z> , 
satisfying the following properties: 

1. the ordering always tells which of two distinct monomials is greater, 

2. the relative order of two monomials does not change when they are each mul- 
tiplied by the same monomial, and 

3. every strictly decreasing sequence of monomials eventually terminates [61]. 
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Different types of monomial orderings exist that satisfy the aforementioned prop- 
erties such as lexicographic, graded lexicographic, graded reverse lexicographic, or 
degree lexicographic. Different monomial orderings are useful depending on the al- 
gorithm that is employed to compute the Grobner Basis. Lexicographic orderings 
will be primarily used throughout our discussion. To learn more about monomial 
ordering, the reader is referred to [52, 60]. 

By fixing a monomial order, we define a leading term for each polynomial of a 
given ideal, denoted as (LT(hi), . . . ,LT(h t )). One can always find a finite subset 
G = (Hi, . . . , H s ) of an ideal I (except for the trivial case / = (0)) such that every 
leading term of / 6 I can be generated by (LT(Hi), . . . , LT(H S )). Here, /el 
means that / is an algebraic combination of hi, . . . , ht, as is required for / to be an 
ideal. Such a subset G is called a Grobner basis with respect to the specific monomial 
order. 3 

One can show that for any given monomial order, every nontrivial ideal / C 
K[xi, . . . ,x m ], has a Grobner basis and that any Grobner basis for an ideal J is a 
basis of /. One can also show that V(I) can be computed by any basis of /, and so 
the solutions of I are the same as that of any of its Grobner basis for any monomial 
ordering. 

A well-defined procedure exists to compute a Grobner basis for any given ideal 
and monomial ordering, called the Buchberger algorithm. It should be noted that the 
Buchberger algorithm reduces to Gaussian elimination in the case of linear equations, 
as it is a generalization of the latter. Similarly, it is a generalization of the Euclidean 
algorithm for the computation of the Greatest Common Divisors of a univariate 
polynomial. 

Comprehensive Grobner Basis 

If the leading coefficient of each element of the basis is 1 and no monomial in any 
element of the basis is in the ideal generated by the leading terms of the other 
elements of the basis, the basis is called a reduced Grobner basis. A reduced Grobner 
basis is unique for a given ideal and monomial ordering, unlike a Grobner basis. 

Now, if we have a parametric ideal, i.e., / = (f\, . . . , f s ) C R[xi, . . . , x n ; ai, . . . , a m ], 
where R is a unique factorization domain, then a Comprehensive Grobner basis 
(CGB) is the distinct reduced Grobner basis for all possible values of the parameters 
ai, . . . ,a m [59]. There are several algorithms available to compute the CGB [62-64]. 
We refer the reader willing to learn more about the actual algorithm and related 
issues to these references, and leave the section by noting that we use the internal 
libraries of Singular to compute the CGB in this paper. 

3 It should be noted here that a Grobner basis may not be unique for a fixed monomial ordering. 
So, we call it a Grobner basis rather than the Grobner basis. However, the so-called reduced 
Grobner basis is unique for a given monomial ordering. The reader is referred to Ref. [52, 60] for 
more details. 
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Let us illustrate with a trivial example of a nonlinear parametric equation ax 2 + 
bx+c = 0; this equation defines an ideal in C[x] with parameters a, b, c in one variable 
x. Singular's CGB library yields that there are three cases: 

1. Case-1: for a = b = c = 0, the solution is the whole C, i.e. all igC; 

2. Case-2: for a = and 6^0, the solution is the line bx + c = 0; and 

3. Case-3: for a ^ 0, the solution is the quadric ax 2 + bx + c = 0, 

as is clearly expected. 

Next, let us consider a more involved example. Take the bi-variate system of two 
equations (fx = ax 2 y 2 + bxy + 2 = 0, f 2 = bx + ay + 2 = 0), where a and b are 
parameters and x and y are variables. Fix the lexicographic ordering x y y. Then, 
the leading terms are clearly LT(f\) = ax 2 y 2 and LT(f 2 ) = bx. The Comprehensive 
Grobner basis is as follows: 

1. Case-1: for a = b = 0, the solution set empty; 

2. Case-2: a ^ 0, b = 0, the solution space is given by ay + 2 = = — 2x 2 — a; 

3. Case-3: a = 0, b ^ 0, the solution space is given by — y + 1 = = bx + 2; 

4. Case-4: ab ^ 0, the solution space is given be —a?y A — 4a 2 y 3 + {b 2 — ^)ay 2 + 
2b 2 y - 2b 2 = and bx + ay + 2 = 0; 

the last three cases can, of course, be checked by simple substitution. 

3 Numerical Algebraic Geometry and Cheater's Homotopy 

Having expounded on the virtues of the CGB, we now introduce a parallel method, 
which attacks our problem from an entirely different perspective. The numerical 
polynomial homotopy continuation (NPHC) method [61] is a recently introduced 
numerical method that finds all the solutions of the given system of polynomial 
equations. It has been used in various problems in particle theory and statistical 
mechanics in Refs. [30-37, 39, 40]. Here we briefly explain the NPHC method to 
make the paper self-contained. 

For a system of polynomial equations, P(x) = 0, where P(x) = . . . ,p m (x)) T 

and x = (xi, . . . , x m ) T , which is known to have isolated solutions, the Classical Bezout 
Theorem asserts that for generic values of coefficients, the maximum number of iso- 
lated solutions in C m is YI4L1 ^ where dj is the degree of the zth polynomial. This 
bound, the classical Bezout bound (CBB), is exact for generic values. The genericity 
is well-defined and the interested reader is referred to Ref. [61, 65] for details. 
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Based on the CBB, a homotopy H(x,t) can be constructed as 

H(x,t) = -y(l-t)Q(x)+tP(x), (3.1) 

where 7 is a generic complex number and Q(x) = (qi(x), . . . , q m (x)) T is a system of 
polynomial equations with the following properties: 

1. the solutions of Q(x) = H(x, 0) = are known or can be easily obtained. Q{x) 
is called the start system and the solutions are called the start solutions; 

2. the number of solutions of Q(x) = H(x, 0) = is equal to the CBB; 

3. the solution set of H(x,t) = for < t < 1 consists of a finite number of 
smooth paths, called homotopy paths, each parameterized by t G [0, 1); and 

4. every isolated solution of H(x, 1) = P(x) = can be reached by some path 
originating at a solution of H(x,0) = Q(x) = 0. 

We can then track all the paths corresponding to each solution of Q(x) = from 
t — to t — 1. The paths which reach P(x) = = H(x,l) are the solutions of 
P{x) = 0. By implementing an efficient path tracker algorithm, all isolated solutions 
of a system of multivariate polynomials system can be obtained: it is shown [61] 
that for a generic 7, there are no singularities (i.e., paths do not cross each other) for 
t G [0, 1). Thus, in the end, we obtain all the solutions of the system P(x) = 0. In 
this respect, the NPHC method has a great advantage over all other known methods 
for finding stationary points. 

Several sophisticated numerical packages well-equipped with path trackers exist, 
such as Bertini [66], PHCpack [67], HOMPACK [68], and HOM4PS2 [69, 70], which are 
all available as freewares from the respective research groups. 

3.1 Cheater's Homotopy 

The advantages of the homotopy based on the CBB are (1) the CBB is easy to 
compute, and (2) the start system based on the CBB can be solved quickly. The 
drawback of it is that the CBB does not take the sparsity of the system into account; 
systems arising in practice have far fewer solutions than the CBB, so a large portion 
of the computational effort is wasted. 

Hence, one can also use homotopies based on tighter upper bounds. For example, 
the 2-Homogeneous Homotopy is constructed by first writing C m = C fc x C k ~ m for 
some k where < k < m, which is accomplished by partitioning the original variables 
into two groups. This has the advantage of incorporating some of the structure of 
the given polynomial system P(x) into the start system Q(x). The corresponding 
bound, called the 2-Homogeneous Bezout Bound (2HomBB), is often tighter than 
the CBB when the polynomial system P(x) has a naturally arising partition of the 
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variables, which occurs in the examples below. Given a partition, the 2HomBB is 
easy to compute and the start system can be solved quickly via linear algebra. 

Another important homotopy is the Polyhedral Homotopy which uses the mono- 
mial structure of the given polynomial system P(x) based on the Bernstein-Khovanskii- 
Kushnirenko (BKK) Theorem [71-73] to yield the BKK bound. Essentially, this up- 
per bound on the number of complex solutions is obtained by computing the mixed 
volume of the convex hull of the Newton polytope (which is based on the exponents of 
the monomials appearing) of each equation. The interested reader from the physics 
community is referred to Ref. [30, 40] for these above two bounds. We note that, 
as with the CBB, the 2HomBB and BKK bounds are also generically sharp with 
respect to the family of polynomial systems under consideration. 

However, in the realistic systems, we do need to take the sparsity of the systems 
fully into account. Indeed, even though the original system may have only a few 
solutions, we may have to track many paths making it computationally expensive. 
The cheater's homotopy is a much more practical method to overcome this diffi- 
culty The crux of cheater's homotopy relies on a theorem, which we state without 
its proof below, that states that the maximum number of solutions of a parametric 
system of polynomial equations over all the parameter-points is the one at a generic 
parameter-point [48]: 

Theorem 1. Let P(x, A) = be a system of polynomial equations, pi(x, A), . . . ,p n {x, A) = 
0, where A = (Ai, . . . , A m ) € C m are parameters and x = (xi, . . . ,x n ) G C n be vari- 
ables. Then, there exists an open, dense, full-measure set U C C n+m such that for 
(6*, ...,&*, A*, ... , X* m ) G U the following holds: 

1. The set X* of solutions x — (xi, . . . , x n ) of the system 



consists ofd isolated points for do < d, where d is the total degree of the system 
for a generic A. 

2. Smoothness and accessibility properties (Properties 3 and 4 of homotopies) still 
hold for the cheater's homotopy which is given as follows: 



Pi(xi, ... ,x n , A*, 



Kn) +K = o, 



Pn\%lj ■ ■ ■ j %>ni ^l; 



\*J + b* n = 



H(x,t) = P{ Xl ,..., 



X n , (1 - t)Xl +*!,...,(!- t)\* m + t\ m ) + (1 - t)b 



where b* = (&*, . . . , &*). It follows that every solution of P(x) = is reached by 
a path beginning at a point of X* . 
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Another way of viewing this is to say that a "special" choice of our coefficients 
may cause the system to be deficient in the maximum number of solutions. If we 
let D be the set of all A that cause these deficiencies, then D is a set of measure 

0. Hence, cheater's homotopy relies on the fact that we choose generic, or random, 
values of the coefficients. If we choose random values for A*, then with probability 

1, A* i D. 

Algorithmically, for a given parametric system P(x, A) = where A = (Ai, . . . , A m ) 
are parameters, we first simply need to solve P(x, A) at a generic parameter-point 
A* G C m . This part has to be solved using a homotopy based on total degree, 
2-homogeneous, or the BKK root counts. Then, in the second part, the system 
P(x, A*) = becomes the start system of all other parameter-points A G C' m — {A*} 
and the solutions of this system become the start solutions. Finally, each path is 
tracked with the below homotopy: 

H(x, A, t) = (1 - t)P(x, A*) + t P(x, A). (3.2) 

The most important trick here is to choose the generic parameter-point from the 
complex space. Once the generic parameter-point is chosen from the complex space, 
the other parameter point can be chosen to be real if the physical situation requires. 
Also note that the gamma trick is implicitly employed since the A* are generic. 

The advantage of the cheater's homotopy over the usual NPHC is huge. In the 
cheater's homotopy, the number of start solutions is drastically smaller than any of 
the other upper bounds usually (in the worst case, it is equal to the smallest of all the 
other upper bounds), and hence the number of paths to be tracked reduces a lot. A 
crucial advantage is that in the cheater's homotopy, we can obtain the start system 
for a generic parametric system once during a computationally expensive 'offline' 
run. This allows us to use the same start system for any number of parameter-points 
in the parameter-space and compute the solutions at each parameter point of interest 
with a much faster 'online' run. 

Moreover, the system at each parameter-point can be solved completely indepen- 
dent from any other parameter-point. In addition, while solving the system at each 
parameter-point, each path can also be tracked independently of all others, making 
it "doubly parallelizable" . We should emphasize that the packages like Bertini and 
PHCpack now have an implementation of cheater's homotopy method in addition to 
the usual NPHC method. 

A recently developed software module of Bertini, called Paramotopy [49], is a dif- 
ferent implementation of cheater's homotopy than the standard one already provided 
with Bertini and can deal with a huge number of parameter points in parallel. In this 
paper we have extensively used this software. Several salient features of this new 
software are: 

1. It offers a few options on choosing the parameter-points: the user can de- 
mand the package to discretize the parameter-space within a given range in as 
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many parameter-points as required or the package can also be fed in a list of 
parameter-points provided by the user. 

2. Paramotopy first stores the data in the RAM before writing to the hard-disk 
which is an efficient data-management practice making the package more effi- 
cient while dealing with many parameter-points simultaneously. 

3. The package is able to disregard any information other than the type of solu- 
tions the user requires in the end, i.e., it can be asked to save only real affine 
solutions, or only non-singular solutions, etc. This is an important aspect of 
the package because by not including all the information produced in a regular 
run, it eliminates massive data proliferation. 

4. If, for any parameter-point where certain paths require higher precision, Paramo- 
topy informs us so that the user can re-run these specific parameter-points at 
higher precision settings. Note that the parameter-points may lie on some al- 
gebraic subset in which exists an algebraic relationship on the parameters. In 
this case, the number of complex solutions shrinks, and hence the related paths 
may require higher precision. 



4 Illustrative Examples 

Having explained our methods in detail, in this section, we will illustrate with ex- 
amples coming from a variety of physical situations, commencing with a toy model 
and moving onto more serious and involved cases. 

4.1 Sysl: A Single-Modulus Example 

We begin with a single-modulus toy example. First, we recall that given the Kahler 
potential K and superpotential W, both as polynomials in fields (pA=i,...,n one can 
proceed to construct the scalar potential from the standard formulae [74]: 



V = e 



K 



K. AB D A WD n W -3IW1 



(4.1] 



As usual the Da represents the Kahler derivative 8a + 9a (K), and K AB is the inverse 
of the field space metric 

JCab = d A d B K . (4.2) 
Our example of this problem requires the solution of the critical set 

d A V = 0, for A = l,...,n . (4.3) 

We can further classify the solutions to (4.3) by the amount of supersymmetry they 
preserve, the value of the bare cosmological constant they dictate, and so forth. The 
most relevant examples are: 
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• SUSY, Minkowski: d A V = 0, D A W = for all A and W = 0; 

• SUSY, AdS: d A V = 0, D A W = for all A but W ^ 0. 

Thus prepared, let us take a single field example, which is addressed in the demo 
of StringVacua and with which we can compare. Let the Kahler potential K and 
superpotential W of an TV" = 1 supersymmetric theory in a single complex moduli 
field T be given as 

K = -3 log(T + f ) , W = a + bT 8 . (4.4) 

Note that the field T comes along with its complex conjugate. Even though they 
can be treated as different variables by merely relabeling them, they are not actually 
independent variables. To avoid this problem, we can write them in terms of real 
and imaginary parts, i.e., T = t + ir where t and r are real. The potential, using 
(4.1), is 

V = ^(46(56(t 2 + r 2 ) 7 - 3a(t 6 - 21tV + 35tV - 7r 6 ))) (4.5) 

Ob 

which has 2 variables. To find the stationary points of V, we need to compute the 
zero locus of the partial derivatives of V with respect to variables t and r: 

^1 = — (46(56(13t 2 - r 2 )(t 2 + r 2 ) 6 - 3a(5t 6 - 63t 4 r 2 + 35t 2 r 4 + 7r 6 ))) = 0, 
at 3t 2 

= l(566r(56(t 2 + r 2 ) 6 + a(9t 4 - 30t 2 r 2 + 9r 4 ))) = . (4.6) 
ot 3t 

For general values of parameters a and b the system already becomes difficult to 
analyze using symbolic methods and one could solve the system for specific values 
such as a = b = 1 [26, 27]. 

We also note that the stationary equations in this example involve denominators. 
Since we are not interested in the solutions for which the denominators are zero, 
we clear them out by multiplying them with the numerators appropriately. In these 
equations, all the denominators are multiples of t. The condition that none of the 
denominators is zero can be imposed algebraically by adding an additional equation, 
1 — yt = 0, with y being an additional variable. Thus, we now have 3 equations in 
3 variables. The CGB library of Singular can deal with this system. The expression 
of the CGB of this system is quite large so we do not write it down here. 

We can also solve this system using Cheater's homotopy. For the range of values 
of a and b given in Figure 1 we find that there are either 6 solutions for a, b > 
and a, b < or 4 solutions for a < 0, b > and a > 0, b < 0; these are indicated 
by different colors in Figure 1. We scan over a total of 100,150 parameter points 
and find a total of 582,676 solutions. Exactly half of the solutions are physical, i.e. 
t > which corresponds to a positive volume of the cycle associated to the Kahler 
modulus T. We did our computation on a desktop machine on single processor (Linux 
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Figure 1. Scanned values for a and b for the one-modulus example in Sysl. The corre- 
sponding number of solutions of the system per parameter point is indicated by the color. 
The total number of parameter points is 100,150 and the spacing between the points is 
equidistant. 

machine with 2.1GHz cloak speed). First, Bertini takes around 30 minutes to solve 
the system from scratch for a given parameter-point. This means that it would have 
taken around 5.731 years to solve the system at all the 100, 150 parameter-points. 
With the CH method, however, we solved them in only 55 hours. 

The scalar masses of the moduli t and r which can be calculated as the eigenvalues 
of the Hessian of V are found to be positive for all parameter points with t > 0. We 
give the moduli masses and the gravitino mass 



which determines the scale of supersymmetry breaking in Figure 2, which is a fre- 
quency plot of the mass values for our space of 100, 150 parameters. 

4.2 Sys2: A Two-Moduli Model 

Now, let us move on to a model with two moduli fields. Consider the Kahler potential 
and superpotential 



ml /2 = e K \W\ 2 , 



(4.7) 



K = -3 log(T + f ) - log(5 + S), 
W = aS + bST + cT 2 , 



(4.8) 
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Figure 2. Masses of the moduli t and r as eigenvalues of the Hessian (left) and the 
gravitino mass (right) for the one-modulus model Sysl. 



with two fields, a Kahler modulus T = t + ir and the axio-dilaton S = s + ia. Hence, 
using (4.1), the stationary equations to be solved are: 

= (3a 2 (s - a)(s + a) + Qat(b(-s 2 + a 2 ) + 2car) + 2bcar(t 2 - 3r 2 ) 

-b 2 (s - a)(s + a)(5t 2 - 3r 2 ) + c 2 (5t 2 - 3r 2 )(t 2 + r 2 )), 
= (-9a 2 (s 2 + a 2 ) + b 2 {s 2 + a 2 )(5t 2 - 9r 2 ) + 2bcr(at 2 + 18str - 9ar 2 ) 

+c 2 (-5t 4 + 2tV - 9r 4 ) + 6a(26(s 2 + a 2 )t + c(hst 2 + Aatr - 3sr 2 ))), 
= (-(cat(6a + bt)) + (6acs + 36 2 (s 2 + a 2 ) - 18bcst - 2c 2 t 2 )r + %car 2 + 6cV), 
= (3a 2 a - 6at{ba + cr) + b(-5bat 2 - ct 2 r + 3bar 2 + 3cr 3 )), 
= 1 - zts . 

The first four equations arise from setting the numerators of the various partial 
derivatives of V to zero, and the last is an auxiliary equation to ensure that the 
denominator does not vanish. Therefore, this system is 5 equations in 5 variables 
with 3 parameters. Again, this is an example used in StringVacua, but a, b and c 
are parameters that were chosen to be 1, —1, 1 respectively in [26, 27]. Now, we can 
undertake the much more challenging task of taking a huge range of parameter- values 
in our computation. 

Again, for generic choice of the values for a, b and c, it becomes difficult for the 
traditional GB method. However, we can compute the CGB for this system; The 
output from Singular of the CGB for this system is quite large, and, hence, we do 
not write it here. 

Using Paramotopy on the other hand, we easily scan over 100, 672 parameter 
points, see Figure 3. We find 6 solutions per parameter point which yields a total of 
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Figure 3. Scanned values for a and b and c for the two-moduli example Sys2. The total 
number of parameter points is 100,672, and the spacing between the points is equidistant. 
The blue regions indicate that there are 6 solutions to Sys2. 

604, 032 solutions. The physicality condition for the case of this model demands again 
t > and also s = g^ 1 > 0, i.e. positive string coupling. This is fulfilled by 503, 299 
solutions. As for Sysl, we did our computation on a desktop machine. Bertini takes 
around 40 minutes to solve the system from scratch for a given parameter-point. 
Hence, for all the 100, 672 parameter-points, it would have taken 7.682 years. How- 
ever, using the CH method, we solved them all in 125 hours. When evaluating the 
Hessian of V, we find that there is always at least one negative eigenvalue, i.e. no 
vacua exist for this model for the set of parameters we considered. 

4.3 Sys3: The Quintic 

An area of string phenomenology where one can make use of the power of Cheater's 
homotopy is the landscape of flux vacua in type IIB string theory [12, 13, 15, 16]. 
As mentioned in Sysl, if one goes to a particular corner of the moduli space, the 
equations 



that determine the supersymmetric vacuum state of the no-scale scalar potential [75, 



are polynomial equations in the complex structure moduli fields 4> a = r, Ui, .., U^.i, 
where r = o + i s (= iS) is the axio-dilaton. The parameters of these equations are 
the flux integers 



DiW = 



(4.9) 



76] 



V = e K K ab D a WD b W , 



(4.10) 
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where F 3 and H 3 are the RR and NS three-form flux of type IIB string theory, and 
(A a , B b ) is a symplectic basis for the b 3 = 2h 2,1 + 2 three-cycles. The flux integers 
are to be chosen freely as long as the D3 tadpole constraint 

I= swi i,lAfl=il,, " y ' (412) 

is not violated. 

For supersymmetric flux configurations, F 3 and if 3 always combine into an imaginary- 
self-dual (ISD) flux G 3 = F 3 — r H 3 [12, 77], which can be written as [9] 

* 6 sH 3 = -(F 3 -aH 3 ). (4.13) 

Consequently, one only has to consider 2h 2,1 + 2 independent directions of the original 
4/i 2 ' 1 + 4 flux integers defining 

*•=(£) md (414) 

Then, the D3 tadpole eq. (4.12) manifests as a positive definite form, i.e. 

L = h 2 + h 2 2 , (4.15) 

for a symplectic basis of three-cycles. With this form of the D3 tadpole, we can use 
Paramotopy to find all flux vacua, i.e. solutions to eq. (4.9) for all flux configurations 
h\, h% with D3 tadpole L < L max for a given maximal tadpole L max . We have to 
take into account the SL(2, Z) invariance of IIB string theory in order to consider 
only physically equivalent flux configurations. It was shown in [78] that for fluxes of 
the form (4.14), only the configurations 

hi \ /— h\\ (— ho\ f ho s 



h 2 J V / \ h J V 

are related by SL(2, Z), and hence physically equivalent. 

This problem to find all supersymmetric flux vacua in the large complex structure 
limit was carried out for the compactification manifold CP 11169 [18] in [78] in the 
context of de Sitter model building in the Kahler uplifting scenario [79, 80]. To 
present a simple example, we consider the one complex structure modulus if) of the 
mirror quintic which has ft, 1 ' 1 = 101 and h 2,1 = 1. In the limit of large complex 
structure, this describes the one moduli subspace of the quintic [81] given by the 
vanishing of the polynomial 

x + x l + x \ + x 3 + x 4 — 5 if) , (4-17) 

in CP* mi with coordinates xq, ..5X4. To obtain a polynomial system let us introduce 
the complex coordinate U = v + iu which relates to if> as 

tf~-Alo g (5V0, (4-18) 

ZTTl 
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up to corrections of 0(U k e 2mU ) that are exponentially suppressed in the large com- 
plex structure limit lm(U) — > oo. This moduli space can be described by an approx- 
imately polynomial prepotential in the large complex structure limit 

5 fcuiV 11 fcoA 2 25wi 25C(3)x' 



with x — 2(/i 1,1 — h 2 ' 1 ) = 200 and u a and Qj, are the periods 



u a = n, g b = n. (4.20) 

J A a JB b 

for a, 6 = 1, 2 and Qb = d^Q. uq can be interpreted as the normalization Q of the 
holomorphic three-form, such that the one physical variable is U = oj\/ujq. Eq. (4.19) 
is valid up to corrections 0(e 2mhU ) that are exponentially small in the limit of the 
large complex structure. We define the large complex structure limit as lm(U) > 2 
such that these corrections to Q are smaller than 10 -3 (?. 4 

Together with the axio-dilaton, r this yields a two moduli example described by 
the Kahler potential and superpotential 

K = K k - log (-i Q(U) A m) ~ (~i(r - f)) , 

= - log '^{QaQa - UJ a Q a )^ - log (-z(t - f )) , 

p / (F 3 — tH 3 ) A Q(U) , 

2 

2VT Yl [(Aa - T h la )G a - (f 2a - T h 2a )U a ] , 



(4.21) 



w 2 



a=l 



where K k is the Kahler potential of the Kahler moduli and setting a' = 1. 

We solve eqs. (4.9) for (4.21) up to a maximal tadpole of L max = 625. This 
corresponds to 481, 825 flux parameter points for which we find a total of 1, 726, 334 
solutions, i.e. on average ~ 3.6 solutions per parameter point. Only 20,280 are 
physical solutions, i.e. gj 1 = s > and u > 2, see Figure 4. If we were to solve it 
at each parameter-point from scratch, then it takes 72 minutes per parameter-point. 
This means that for total 481, 825 parameter-points that we solved the system for, we 
would have taken 481, 825 x 72/60 = 578, 190 hours. But instead, using Paramotopy, 
we solved all of them in only 3776 hours, i.e., around 26.4 human hours with a cluster 
of 144 processors. We could not obtain the CGB for this system in a reasonable time. 



4 Note that there is a conifold singularity at V> 5 = 1 that is excluded from the large complex 
structure limit Im(U) > 2 according to eq. (4.18). 
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Figure 4. Left: Distribution of all algebraic solutions to eq. (4.9) with s > and u > 
for the Kahler potential and superpotential (4.21) in Sys3. The gray shaded area indicates 
the space of physical solutions, i.e. s > and u > 2. Middle: Distribution of the dilaton 
r = a + i s. Right: Distribution of the superpotential Wo- 

We can make use of the SX(2,Z) symmetry of IIB string theory, to transform 
each solution to the fundamental domain 

- \ < Re(r) < \ and |r| > 1 , (4.22) 

via the successive transformations 

t' = t + &, W^ = W , and t' = -l/r, = W /r (4.23) 

with b G Z. We show the distribution of the obtained values for r = a + is and 
the flux superpotential Wo in Figure 4. We see that the the strongly coupled region 
s = l/g s ~ 1 is preferred and values of Wo ~ C(10 2 — 10 3 ) are preferred. The same 
qualitative behavior was observed in [78] for the manifold CP 11169 [18]. 

Finally, we give the masses m 2 of the moduli t, r, s and a and the gravitino mass 
mL in Figure 5. For this, we have to specify the value of the Kahler moduli 
Kahler potential in eq. (4.21). For a Calabi-Yau compactification, this is given as 
Kk = —2 log V, where V is the volume modulus of the Calabi-Yau and its vacuum ex- 
pectation value depends on the stabilization mechanism of the Kahler moduli. Here 
we choose V = 100 in string units for definiteness. 

4.4 Further Applications: Sys4 

As mentioned in the introduction, large systems with a multitude of parameters are 
ubiquitous. Thus, our methods above should be applicable to far more situations than 
finding the extrema of a PEL in the context of effective supersymmetric Lagrangians. 
In this subsection, as a parting example, let us see the power of Cheater's homotopy 
applied to the geometry of Calabi-Yau manifolds. 
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Figure 5. For Sys3 ((4.9) and (4.21)), the distribution of the masses m 2 of the moduli 
t,r,s and a, i.e. all eigenvalues of the Hessian of V (left) and the gravitino mass m^ 2 
(right) in units of Mj for V = 100. 



For concreteness and continuing along the vein of Sys3, let us focus on the problem 
of finding the singular locus (or the absence thereof) given a family of quintics. 5 

Let us consider the family of quintic manifolds given by as a homogeneous hyper- 
surface in CP 4 with coordinates xq, . . . , X4, and let 

Q(xo, . . . , X4; a, b, c) = Xq + x 1 + x\ + x\ + x\ + 0x0X1X2X3X4 + 6x 4 xi + CX1X2X3, (4.24) 

where a, b, c are complex parameters. Of course, the most general family of quintics 
has 101 deformation parameters corresponding to h 2,1 of the manifold, but this ex- 
ample suffices to show the power of the method. Now, for b = c = 0, we have the 
one most familiar to us, with a = —5 being the conifold point in the complex struc- 
ture moduli space. The singular locus of Q is given by the Jacobian ideal (d Xi Q) 
for i = 0,...,4, excluding the origin and defined up to scaling of the projective 
coordinates. 

It is more convenient to work over affine patches of CP 4 which can be rescaled to 
be Pi = {xi = 1} for % = 0, . . . , 4. For each of the five patches, we compute the four 
remaining partial derivatives, the solution of which is then the system we need to 
analyze. That is, we need to perform Paramotopy on 

d xm Q{xi = 1; a, b, c) = , z = 0,...,4. (4.25) 

5 Recently there has been much activity in study the database of so-called complete intersection 
Calabi-Yau (CICY) manifolds of which the quintic is the simplest example as well their smooth 
quotients. Checking smoothness for models with this database of 7890 Calabi-Yau threefolds and 
their descendents, for example, is a crucial step [82-84]. 
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Figure 6. The number of solutions at each parameter point {a, b, c} for system Pq (top 
left), Pi (top right), P2 and P3 (bottom left) and P4 (bottom right) indicated by the 
different colors. 

All of the 5 systems are of the same size and we take a, b, c G [—25, 25] in increments 
of one, i.e. 51 3 = 132,651 parameter points. Bertini takes from 144 to 223 seconds 
at a parameter-point if solved from scratch, hence it would have taken from 0.6057 
to 0.938854 years for all 132651 parameter-points. Paramotopy solved each of these 
systems for all parameter-points in from 19.216 to 92.9002 hours (using 72 processors, 
it took 960.812 to 4645.01 seconds). Obtaining the CGB for this system is quite fast 
as well; however, obtaining the solutions for each parameter-point from the CGB 
takes the same amount of time as Paramotopy. Hence, we show the results only 
obtained from Paramotopy here. We show the number of solutions per parameter 
point indicated by the different colors in Figure 6 and Table 1. The solution space 
of the systems P2 and P3 is identical since eq. (4.24) is invariant under x 2 ^ ^3- 

5 Conclusion and Outlook 

Parametric systems of non-linear equations arise very naturally and frequently in 
theoretical physics. The parameters add one more hurdle in solving such systems in 
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#soZ 





1 - 10 


11 - 100 


101 - 240 


Po 


35,613 


93, 162 


1,250 


2,626 


Pi 





109,326 


20, 724 


2,601 


P2,Pl 





126, 159 


6,441 


51 


Pa 


60 


129,990 





2,601 



Table 1. Distributions of the number of solutions j^sol for the five systems P{. 

addition to the non-linearity of the equations as one has to solve the system many 
times. In the present paper, we have demonstrated that two of the sophisticated 
algebraic geometry methods can solve such systems very efficiently provided that 
the nonlinearity is polynomial-like. 

The first method is the Comprehensive Grobner basis (CGB) which is a symbolic 
method based on Grobner bases. Given a (non-parametric) system of equations, 
we can systematically find another system of equations which has the same solution 
space as the original one and is easier to solve, called a Grobner basis. If we have a 
parametric system of equations, then we can compute the CGB which is a Grobner 
basis for all values of the parameters including their special values. To obtain the 
solutions of the system at a specific parameter-point, one then needs to input the 
values of the parameters in the CGB and solve the system using the traditional 
methods. However, the algorithmic complexity to compute the CGB is the same, or 
in many cases even worse, as the Grobner basis method, i.e., the memory required by 
the algorithm blows up exponentially with increasing number of variables, equations, 
monomials, and degree. Moreover, the related algorithm is highly sequential and 
unparallelizable. Hence, though the CGB method is extremely useful in the cases 
where we can finish the computation, but for many physical systems the method 
often falls short. 

The second method we described in this paper is called the cheater's homo- 
topy. Here, one first solves the given parametric system of polynomial equations 
at a generic parameter-point because it is shown that the number of solutions at 
such a parameter-point is an upper bound on the number of solutions at any other 
parameter-points. Then, one tracks solution-paths from the solutions at the generic 
parameter-point to the ones at which the system needs to be solved using the numeri- 
cal polynomial homotopy continuation method. This method is highly parallelizable, 
and, hence, we can solve the system at thousands of parameter-points in a short time 
on a computer cluster. We use a package called Paramotopy which is not only a highly 
efficient implementation of the related algorithm, but also does data and memory 
management effectively. These features allow for a huge number of parameter-points 
to be solved simultaneously. 

We have used these two methods on various examples arising from computing the 
extrema of potential energy landscapes, especially in the context of moduli stabiliza- 
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tion in string phenomenology. However, obtaining the solutions for each parameter- 
point from the CGB (for the cases it is possible to obtain the CGB) takes the same 
amount of time as Paramotopy; hence, we show the results only obtained from the 
cheater's homotopy in this paper. 

The examples are given as a 4D M = 1 effective supergravity scalar potential, 
defined by a Kahler potential K and a superpotential W . Our first two examples 
Sysl and Sys2 are toy models of one and two complex scalar fields, respectively. 

We looked for stationary points of these scalar potentials for ~ 10 5 parameter 
points, finding half of the stationary points are vacua in the case of Sysl and finding 
no vacua in the case of Sys2, i.e., at least one tachyonic direction always exists. 
Finally, we looked at flux compactifications on a realistic Calabi-Yau, the quintic in 
the case of two moduli. We look for flux vacua for ~ 5 ■ 10 5 flux configurations, which 
corresponds to a maximal D3 tadpole of L = 625. Our results indicate preferences of 
strongly coupled vacua g s > 1 and values of the superpotential Wo ~ C(10 2 — 10 3 ). 

The applicability of the methods described here is, of course, much wider and is in- 
tended for any large parametric polynomial systems, with the parameters appearing 
explicitly as coefficients. Indeed, these methods should provide an extremely pow- 
erful tool for many areas in theoretical physics such as potential energy landscape, 
statistical mechanics, particle phenomenology, string phenomenology, and non-linear 
dynamics. As an illustration, we have applied Cheater's homotopy to the geometry 
of Calabi-Yau manifolds in the case of Sys4- 

If we were to solve Sysl, Sys2, the quintic, and the Calabi-Yau systems from 
scratch for each parameter-point using the numerical homotopy continuation method, 
a single processor machine would have taken 5.731 years, 7.682 years, 66.003 years, 
and 0.6057-0.938854 years, respectively. But using the cheater's homotopy method 
we solved them in 55 hours, 125 hours, 3776 hours and 19.216-92.9002 hours, respec- 
tively. This is a drastic reduction of the computational efforts, see also Table 2. 



System 


Number of parameter points 


single-processor 


Cheater's homotopy 


Sysl 


100, 150 


5.731 yr 


55 h 


Sys2 


100,672 


7.682 yr 


125 h 


Sys3 


481,825 


66.003 yr 


3776 h 


Sys4 


132,651 


0.6 - 0.9 yr 


19.2 - 92.9 h 



Table 2. Number of parameter points and computation time for the examples under 
consideration. 



We emphasize that the two methods are based on complex algebraic geometry 
which means that the variables are first considered to be complex, and then only 
purely real solutions are retained, if desired. Because in many applications we only 
have real parameters and only real solutions are physically important, a natural ques- 
tion to ask is if any method exists which can deal with systems with both parameters 
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and variables defined over reals. Indeed, a more recent method exists from symbolic 
real algebraic geometry, based on the discriminant varieties [50, 51], which treats 
both variables and parameters as reals from the beginning and gives a different set 
of information than the methods described in this paper. This certainly constitutes 
a direction worthy of pursuit. On the numerical side, little progress exists in this 
direction, and so far methods exist to extract only limited information out for the 
whole parameter space. For example, in Ref. [85], a numerical method is proposed 
which can tell us the maximum and minimum number of real solutions over the 
whole parameter-space. It is our hope that the applications in this paper will also 
become a motivation to develop such methods and making them more for real-life 
applications. 
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